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We study the properties of normal, superconducting (SC), and CDW states for an attractive Hubbard model on the 
square lattice, using a variational Monte Carlo method. In trial wave functions, we introduce an interspinon binding 
factor, indispensable for inducing a spin-gap transition in the normal state, in addition to the onsite attractive and intersite 
repulsive factors. It is found that, in the normal state, as the interaction strength \U\/t increases, a first-order spin-gap 
transition arises at |t/cl ~ W (W: bandwidth) from a Fermi liquid to a spin-gapped state, which is conductive as a 
result of the hopping of doublons. In the SC state, we confirm by the analysis of various quantities that the mechanism 
of superconductivity undergoes a smooth crossover at approximately \Uco\ ^ Wd from a BCS type to a Bose-Einstein 
condensation (EEC) type, as \U\/t increases. For \U\ < \Ucol quantities such as the condensation energy, a SC correlation 
function and the condensate fraction of onsite pairs exhibit the behavior of ~ Qxp(-t/\U\), as expected from the BCS 
theory. For \U\ > |t/coL quantities such as the energy gain in the SC transition and superfluid stiff'ness, which is related to 
the cost of phase coherence, behave as ~ ^^/|f/| oc Tc, as expected in a bosonic scheme. In this regime, SC transition is 
induced by a gain in kinetic energy, in contrast to the BCS theory. We refer to the relevance to the pseudogap in cuprate 
superconductors. 

KEYWORDS: superconductivity, BCS-BEC crossover, pseudogap, superfluid stiffness, condensate fraction, 
CDW, attractive Hubbard model, square lattice, variational Monte Carlo method 



1. Introduction 

Since the pioneering effort by Eagles,^ many researchers 
have extensively and repeatedly addressed the transition of 
superconducting (SC) properties from a BCS type to a Bose- 
Einstein condensation (BEC) type as the strength of attrac- 
tive potential between fermions is increased. Following early 
studies of the BCS-BEC crossover^' ^ using continuum mod- 
els with the superfluidity of ^He in mind, Nozieres and 
Schmitt-Rink^ showed by approximation that the SC prop- 
erties smoothly evolve with the correlation strength in an at- 
tractive Hubbard model (AHM). Later, stimulated by the dis- 
covery of high-Tc cuprates with a small coherence length, 
numerous researchers have tackled AHM,^ especially in two 
dimensions (2D); now, in connection with the evolution of 
pseudogaps as the doping rate S decreases in the so-called 
underdoped regime,^ the problem of BCS-BEC crossover as 
a function of ^ is a subject of urgency.^' ^ Entering this cen- 
tury, we have become capable of directly observing phenom- 
ena of crossover^' ^^ and pseudogaps^ ^' ^^ in traps of ultracold 
dilute alkali gases, ^^' ^^ for which physical parameters can be 
artificially tuned. Recent experimental advances have brought 
hope of obtaining similar observations on optical lattices. 

In the above stream of research, AHM is one of the most 
important and basic lattice models for studying the evolu- 
tion of SC properties according to the interaction strength 
U/t (U: onsite interaction strength, t\ hopping integral be- 
tween nearest-neighbor sites). In early and later studies of 
AHM, mean-field-type^ and diagrammatic^' ^^~^^ approaches 
were used; although they successfully treated the weakly 
correlated regime, where the original BCS theory is basi- 
cally valid, and developed a conceptual framework of the 
BCS-BEC crossover, more reliable methods remain neces- 
sary to establish properties in the intermediately (unitary) and 
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Strongly correlated (BEC) regimes. First, as an unbiased way, 
quantum Monte Carlo (QMC) calculations were implemented 
in the weakly and intermediately correlated regimes (| [/| < W, 
W: bandwidth), because QMC is free from the negative-sign 
problem for AHM, but statistical fluctuation increases with 
increasing \U\/t and system size. In 2D, SC transition of the 
Berezinskii-Kosterlitz-Thouless type was confirmed and an 
^-dependent phase diagram was discussed (n\ particle den- 
sity). ^^ Then, it was shown in the normal state (T > Tc) for in- 
termediate \U\lf^ that a thermal-activation-type behavior ap- 
pears in magnetic susceptibility, but that the charge compress- 
ibility is almost T-independent.^^ This spin-gap behavior was 
corroborated by a peak split in the density of state. ^^ The dy- 
namical mean field theory (DMFT), which becomes exact in 
infinite dimensions and is applicable to an arbitrary interac- 
tion strength, is another important approach to AHM. Early 
DMFT studies addressed normal branches without introduc- 
ing SC orders at low temperatures {T < T^^^ and T = 0^^), 
and found that the normal state undergoes a first-order transi- 
tion from a Fermi liquid to a gapped state at|t/|/W= 1-1.5, 
as \U\/t increases. Later, various properties of the SC phase 
were calculated,^^"^^ and the crossover was characterized by 
the SC gap and superfluid stifl'ness.^^ 

Another efl'ective approach to AHM is a many-body vari- 
ation theory, which is applicable continuously in the entire 
range of correlation strengths and particle densities. In con- 
trast to DMFT, the dimension and lattice form are realistically 
specified, and one can treat wave-number-dependent proper- 
ties in low-lying states. Furthermore, since wave functions 
are explicitly given, this approach has advantages in forming 
a physical picture. Because an AHM of a bipartite lattice is 
mapped to a repulsive Hubbard model (RHM) by a canonical 
transformation,^^"^^ one can develop a theory relying on the 
knowledge of RHM. Thus, the well-known Gutz wilier wave 
function (GWF)^^ became a primary trial function for the 



J. Phys. Soc. Jpn. 



DRAFT 



normal state; first, its properties were studied ^^ using the so- 
called Gutzwiller approximation (GA).^^ As known for RHM, 
although GWF itself is always metallic, ^^ additional GA in- 
duces a spurious metal-insulator (Brinkman-Rice) transition^"^ 
at finite I^/brI/^ in finite lattice dimensions. For \U\ > I^/brI 
in AHM, all the particles tightly form onsite singlet pairs, 
and hopping completely ceases, so that the Brinkman-Rice 
transition remains a metal-insulator (Mott) transition also in 
AHM. Later, approximations similar to GA, which may be 
correct in infinite dimensions, have also been applied to the 
SC state^^"^^ to discuss the BCS-BEC crossover. However, 
to avoid the ambiguity of GA in realistic dimensions and to 
make use of the merits of the variation method, we need to ac- 
curately estimate variational expectation values. This claim is 
satisfied by a variational Monte Carlo (VMC) method,^^'^^^^ 
which treats local correlation factors exactly without a minus 
sign problem. A decade ago, a VMC method was applied to 
a normal state in AHM to study a transition corresponding 
to the Mott transition in RHM by introducing a binding fac- 
tor between adjacent antiparallel spinons."^^ For simplicity, we 
call a singly occupied site a spinon. However, the interpreta- 
tion of the transition was incorrect on account of the limitation 
of treated system sizes and an insufficient analysis. Recently, 
VMC has been applied to solving problems with optical lat- 
tices in a confinement potential. "^^ 

In this study, on the basis of VMC calculations of high pre- 
cision for normal, SC and CDW states, we modify the pre- 
vious results"^^ and make features of the BCS-BEC crossover 
in AHM on the square lattice microscopically more clear. We 
mainly discuss the following points: (1) In both normal and 
SC states, a correlation between adjacent antiparallel spinons, 
in addition to the Gutzwiller correlation, is indispensable to 
qualitatively derive proper behavior. (2) In the normal state, 
which underlies the SC state, a first-order phase transition oc- 
curs at |t/cl ~ W from a Fermi-liquid to a spin-gapped state. 
This transition is caused by the competition between the size 
of an antiparallel- spinon pair and the interpair distance, as in 
the case of Mott transitions in RHM.^^'^^ (3) The properties 
of SC noticeably change at approximately \Uco\ ~ \Uc\, which 
are compared with those derived in a strongly correlated RHM 
for high-Jc cuprates. Part of the present result has been pub- 
lished before.^^ 

The rest of this paper is organized as follows: In §2, we 
explain the model and method used in this study. In §3, we 
provide a discussion of the spin-gap transition arising in the 
normal state, and of the features in the spin-gapped regime. In 
§4, we consider a BCS-BEC crossover from various points of 
view. In §5, we briefly summarize our main results. 

2. Formulation 

In §2.1, we introduce AHM and briefly mention its relation 
to RHM. In §2.2, we discuss trial wave functions for normal, 
SC, and CDW phases. In §2.3, we briefly explain the setup of 
VMC calculations in this work. 

2. 1 Attractive Hubbard model 

We consider a single-band attractive Hubbard model {U < 
0) on a square lattice: 



where ni 



ko- / 



^jT^iU' 



(1) 



ijcr - c^jo-^jo-, Cjo-, and c^o- are fermion annihilation 
operators in the Wannier and Bloch representations, respec- 
tively, and 



si^ = -2t(coskx + cos^y). 



(2) 



We use the hopping integral t and lattice constant as the units 
of energy and length, respectively. Because the lattice has a 
particle-hole symmetry at n = N/N^ = 1 (N: number of par- 
ticles, A^s- number of lattice sites), and properties at half fill- 
ing are deduced from the results of RHM using correspond- 
ing wave f unctions, ^^'^^ as mentioned below, we mostly treat 
cases of ^ < 1. The chemical potential term -^Yjja^jcr may 
be added to eq. (1) to adjust particle density, if necessary. 

In the following, we summarize the relation of AHM to 
RHM. The attractive Hubbard Hamiltonian eq. (1) on a bi- 
partite lattice satisfying the relation ^k = -^-k+Q [here Q = 
(n, 7i)] is mapped by the canonical transformation^^' ^^ 



<^kT = <^kT' <^ki 
to RHM with constant shifts: 



■k+Qi 



(3) 



^=Z^^^kcr^^^ + 



ko- 






where hjrr = cl^Cja-, S^ - 



^4) 
jcr - ^jcr^ja, ^ j - {rij^ - hji)/2 and /z^ = Ncr/N. 
A tilde denotes the representation transformed according to 
eq. (3). The chemical potential f and n in AHM are related to 
the eff'ective magnetic field as /z = 2^ and to the magnetization 
as m = 1 - ^ in the z-direction in RHM, respectively. There- 
fore, unless the original AHM has a spin polarization (m = 0), 
the particle density in the transformed RHM is always at half 
filling (h = 1). Also, the order parameters of CDW and onsite 
singlet pairing defined as 



^CDW 



1 

N 



Y^e^Q-^^{nj,^nj,-l} 



(5) 



in AHM are transformed into the forms of the z- and xy- 
components of the SDW order parameter: 



Oi 



SDW 



n- 



1 

N 



Y^e^^-^^{nj^-njO 



l^<cj^c,,) or -;^( 



ji^J 



(7) 



c,T>, (8) 



respectively, in RHM.^^ It is widely accepted that, at T = 0, an 
antiferromagnetic (AFM) long-range order with equal magni- 
tudes of O^j^^ for a = x,y,z arises in the half-filled RHM on 
the square lattice for arbitrary U (> 0), and that the AFM or- 
der in the z-direction is easily destroyed by a field applied in 
the z-direction h, whereas the AFM orders in the xj-plane sur- 
vive. This implies that the ground state of AHM possesses a 
singlet pairing order for any U and n (f ), and simultaneously 
possesses a CDW order of the same magnitude at half fill- 
ing n = I (^ = 0)?^ This argument was confirmed by direct 
calculations for AHM.^ Although the above mapping holds 
unconditionally in exact treatments, when some approxima- 
tion is applied, the validity of the mapping has to be verified 
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individually for each specific treatment. 

2.2 Trial wave functions 

As a development of our previous study, ^^ we apply a 
many-body variation theory to the Hamiltonian eq. (1). As 
a trial wave function, a two-body Jastrow-type ^ = ^Omf 
was adopped,"^^ where Omf is a one-body (mean-field) wave 
function and ^ is a many-body correlation (Jastrow) factor. 

As the many-body part, we use the form V - VfPqPQ in 
this work. The onsite (Gutz wilier) projector 



Vo^W\\-{\-g)d\^ 



(9) 



(10) 



with dj - nj^nji, is the most important. The variational pa- 
rameter g increases the number of doubly occupied sites (dou- 
blons), and ranges over 1 < g < oo for U < 0; in the limit of 
g ^ oo, singly occupied sites (spinons) are not allowed in a 
nonmagnetic case. If we put g = l/g, the properties of Pcig) 
for RHM are applicable to the present case Pcig)-"^^ 

To explain the importance of a binding factor between the 
up and down spinons Pq, it is convenient to refer to an eff'ec- 
tive Hamiltonian in the strong-correlation limit (t/\U\ -^ 0):^^ 

with 

bi = Ci^Cii, Pi = -(m^ -h mi - 1), (Ti = -(m^ - nn). (11) 

The first term of eq. (10) indicates the hopping of doublons. 
The second is a repulsive interaction between doublons [or 
empty sites (holons)] and an attractive interaction between a 
doublon and a holon in nearest-neighbor (NN) sites. The third 
works as an AFM-Ising interaction. The expectation values of 
these terms can be reduced using antiparallel-spinon config- 
urations in NN sites. To encourage such configurations, we 
introduce the attractive intersite correlation,"^^ 



J 



/ ) 



(12) 



(13) 



where s^ = njo-(l - nj-o-) (spinon projector), and r runs over 
NN sites of the site j. In Pq, the parameter yu (0 < yu < 1) con- 
trols the strength of binding between NN antiparallel spinons; 
for yu = 0, spinons are free of binding, while in the limit 
yu ^ 1, antiparallel spinons are necessarily paired as near- 
est neighbors. As we will see later, Pq is indispensable for a 
spin-gap transition"^^ and a proper description of the SC state. 
In fact, Pq is the canonical transformation through eq. (3) 
of the doublon-holon binding factor often used to describe 
Mott transitions in RHM.^^'^^ A Mott transition in RHM cor- 
responds to a spin-gap transition in AHM, as we will see in 
§3. Since Qj is a spin-dependent projector, the so-called spin 
contamination^^ arises in the wave function, namely, ^ de- 
viates from an eigenstate of S^ = (2j Sy)^. However, in this 
case, the expectation values of (S^) estimated using a VMC 
method are as small as 0.15 (2) in the SC (normal) state for 
N = 200-300 at its maximum at U - W. Because these val- 
ues, particularly of the SC state, are much smaller than those 
of the AFM state, the spin contamination is considered to have 



little influence on the results. 

As a factor supplementary to Pq, sl repulsive correlation 
suited to eq. (10) should be considered. As a simple one, we 
check a repulsive factor between NN doublons: 



Pf 



= Y\\i-fdj{i-Y\djJ 

/ L V T / 



(14) 



where /(0</< l)isa parameter, dj = 1 - dj, and r 
runs over NN sites of the site j. The projector Pf reduces the 
weight of configurations with adjacent doublons by 1 - /; for 
/ ^ 0, the eff'ect of Pf vanishes, and for / ^ 1, a doublon 
cannot sit in a NN site of another doublon. 

Now, we turn to the one-body part O of the wave function. 
For a normal state, we adopt the Fermi sea Op. Since general 
features of ^n = P^f with / = were studied in a previ- 
ous paper, "^^ here we focus on the properties of the transition 
arising atU-W, which was regarded as a Mott transition. "^^ 

It is known that the BCS state Obcs can deal with the BCS- 
BEC crossover in some degree;^' ^ it is natural to employ Obcs 
for a SC state: 



o 



BCS 



= ^ akcl^c"^ 



kT -ki 



10). 



where the particle number is fixed and 

Ap 



Vk 

ak = — 
Wk 



^k - ^ + -^(6:k - 0^ + A^ 



(15) 



(16) 



with Ap and ^ being variational parameters corresponding to 
the SC gap Asc and chemical potential ^, respectively, in the 
weakly correlated limit, and 



"k (v^) = i 



1 + (-). 



Sk-C 



^I{^^ 



+ A^ 



For Ap -^ 0, Obcs is reduced to Op. Here, we assume Ap to 
be a homogeneous s wave on account of the attractive contact 
potential. For RHM, a form similar to eq. (16) with a dx2_y2- 
wave pair potential was studied,"^^' ^^ where, as S decreases, 
what Ap means deviates from the SC gap, and represents a 
pseudogap.^"^'^^ In contrast, in the present case, Ap seems to 
reflect the magnitude of T^ for any \U\/t, except for ^ -^ 1, 
for which T^ is considered to vanish owing to the CDW order. 
The correlated SC function ^sc = ^Obcs is mapped using 
eq. (3) to a projected AFM wave function ordered in the x-y 
plane at ^ = 1 . 

In addition, we check a CDW wave function for ^ ~ 1 : 

OcDw = n (-^kcL +^k4Q.) I0>' (17) 

k,o- 

where the k sum is taken in the Fermi sea, Q = (n, n), and 



c^k (13k) = 



A 



!-(+)■ 



£k 



^ 



^k+^? 



(18) 



with Ac being a parameter corresponding to the CDW gap. 
^CDW = ^OcDW is mapped through eq. (3) to a projected 
AFM wave function ordered in the z-direction at ^ = 1. Be- 
cause the AFM order is isotropic in RHM, ^sc and ^cdw 
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should yield identical results at half filling. 

2.3 Variational Monte Carlo method 

In estimating variational expectation values with respect to 
^ discussed in §2.2, we use a VMC method,^^'^^^^ which 
gives virtually exact values for finite but relatively large sys- 
tems. Since the number of variational parameters is not large, 
we execute rounds of linear optimization for each parameter 
with the other parameters fixed until the parameters as well 
as the energy converge (typically 3-5 rounds) with 2.5 x 10^ 
particle configurations generated through a Metropolis algo- 
rithm. After the convergence, we continue to execute addi- 
tional 20 to 30 rounds of iteration with successively renewed 
configuration sets. We determine the optimized values by av- 
eraging the data obtained in the additional rounds; in aver- 
aging, we exclude scattered data beyond the range of twice 
the standard deviation. Thus, the optimal value is an average 
of substantially more than several million samples. Physical 
quantities are computed with the optimized parameters thus 
obtained with 2.5 x 10^ samples. 

We use systems of LxL sites of up to L = 32 for ^n and L - 
24 for ^sc with the periodic-antiperiodic boundary conditions 
to reduce level degeneracy. We choose the particle densities to 
satisfy the closed-shell condition, and mainly study n - 0.25 
(0.26), 0.5, and 0.75. 

3. Spin-Gap Transition in Normal State 

As mentioned in §2.1, the ground state of AHM is SC for 
any Ujt and n (and CDW at ^ - 1). Therefore, the normal 
state, we address in this section, is not the ground state of 
eq. (1). The significance to study ^n is not only in a passive 
sense that a normal state appears when the SC state is de- 
stroyed by, e.g., magnetic field or impurities, but in that ^n 
underlies ^sc, just as a SC transition is understood by the 
instability of the Fermi sphere against an infinitesimal attrac- 
tive interaction in the BCS theory. Namely, normal states are 
deeply involved in the mechanism of SC transitions. 

In §3.1, we show the improvement in energy by introduc- 
ing Pq and ^/, and the energy gains using the SC and CDW 
states. In §3.2, we confirm the existence of a spin-gap transi- 
tion in the normal state. In §3.3, we consider the mechanism 
of the spin-gap transition and other properties. 

3. 1 Energy improvement 

First, we briefly look at the energy improvement by the pro- 
jection factors Pq and Vf in the normal, SC, and CDW states. 



Table I. Comparison of energy per site Elt among three kinds of projec- 
tion factors pQ, 9q = PqPq, and Py^ = P/PqPq in normal and SC states. 
A system ofn = 0.5 and L = 12 is used; the tendency is the same for other 
n's. The last digit in each section includes some statistical errors. 



\u\/t 


1 


3 


7 


10 


Pg^f 


-1.37513 


-1.54791 


-2.0677 


-2.6215 


PQ^F 


-1.37536 


-1.55099 


-2.10435 


-2.7417 


Pr^F 


-1.37536 


-1.55097 


-2.10440 


-2.7416 


Pg Obcs 


-1.37531 


-1.55686 


-2.17815 


-2.80500 


"Pq <I>bcs 


-1.375535 


-1.55908 


-2.18444 


-2.81162 


"Pr ^Bcs 


-1.375534 


-1.55918 


-2.18681 


-2.81436 



normal state is considerably improved by Pq on that of GWF, 
especially for large \U\/t's (Table I). Moreover, it is known 
that the phase transition discussed in §3.2 does not arise with- 
out Pq.^^'^^ Thus, the factor Pq is indispensable to appro- 
priately describe the normal state. These aspects of Pq cor- 
respond to those of the doublon-holon (D-H) binding factor 
in RHM."^^'^^ On the other hand, the improvement by Pf on 
PqPg^ is almost imperceptible for any \U\/t as shown in Ta- 
ble I. The optimized parameter / is nearly zero, namely, Pf 
scarcely modifies the wave function. 

In the SC state, the improvement in E/t by Pq on ^g^bcs 
is not as large as that for the normal state. This is because the 
eff'ect of binding between up and down spinons, i.e., the eff'ect 
of singlet pair creation, is already included in Obcs to some 
extent, as the one-body AFM state has some D-H binding ef- 
fect for RHM.^^ Further energy reduction by Pf is again neg- 
ligible for small \U\/t's and remains relatively small in magni- 
tude for larger \U\/t's, as compared with the energy reduction 
by Pq. The magnitude of energy reduction by Pf is similarly 
smaU for ^ = 0.25 and 0.75. 

Since we find that a short-range repulsive factor Pf pro- 
duces only negligible eff'ects in all the cases we treat, we omit 
Pf and use the form P = Pq = PqPq as the many-body fac- 
tor in ^ in the rest of this paper, unless otherwise specified. 



L=12 



Wi ^SC -^CD W 




Fig. 1. (Color online) Particle-density dependences of the energy gains by 
^sc and ^cdw are plotted for five values of \U\/t. For small \U\/fs, the 
system size dependence (fluctuation) is large owing to the large coherence 
length. 

Finally, we compare the energy gains using the SC and 
CDW states: 



A£'sc (A£'cDw) - E^ - Esc (£'cdw), 



(19) 



where £'n, E^c, and Eqdw are the optimized energies per site 
for ^N, ^sc and ^cdw, respectively. Figure 1 shows the n de- 
pendences of A£'sc and A^cdw for large ^'s. At half filling 
(n = 1), the SC and CDW states are degenerate, but this de- 
generacy is immediately lifted for ^ < 1 . A£'cdw rapidly dete- 
riorates and vanishes as n decreases, whereas A^sc preserves 
appreciable values for high densities and gradually decays un- 
til ^ = (not shown). This feature of AE coincides with what 
we discussed for the canonical transformation in §2.1.^ 
In the remainder of this section, we will concentrate on ^n- 



As studied in detail in ref. 41, the variational energy in the 
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3.2 Spin- gap transition 

In a previous VMC study using Pq^j^,"^^ a transition was 
detected with systems up to L = 12 at U = U^ - 9t. How- 
ever, this transition was misinterpreted as a continuous metal- 
insulator transition. In this subsection and the next, we study 
the features of this transition more carefully. 



normal 




r(0,0) 



X(7r,0) M(7i,7r) 

k 



r(o,o) 



Fig. 2. (Color online) Momentum distribution function for some values of 
\U\/t near transition point \Uc\/t along path (0, 0)-(7r, 0)-(;r, ;r)-(0, 0). The dis- 
continuity at kp indicated by an arrow is used to estimate Z shown in Fig. 3. 

First, we confirm the existence of a transition. In Fig. 2, we 
plot the momentum distribution function 



""^^^ = 9 Z (^cr^kcr) 



(20) 



for \U\/t ~ 9 and ^ = 0.25 (L = 32). For \U\/t < 9.0, ^(k) 
has discontinuities on the F-X and F-M segments, indicating 
that a Fermi surface exists and the state is a Fermi liquid. 
On the other hand, the discontinuity suddenly vanishes for 
\U\/t > 9.05, and ^(k) becomes a smooth function of k. It fol- 
lows a certain gap opens and the state becomes a non-Fermi 
liquid for |^| > \Ucl with 9.0 < \Uc\/t < 9.05 in this case. 
Through similar analyses, we found \Uc\/t ~ 0.875 (0.83) 
for n = 0.195 (0.121) for L = 32; \Uc\/t tends to gradually 
decrease with n.^^ Thus, a transition from a Fermi liquid to 
a non-Fermi liquid certainly exists, as found in our previous 
study.^^ According to similar analyses for L = 24 and 28 and 
n ~ 0.25, the system- size dependence of \Uc\/t is very small 
at these values of L, but \Uc\/t tends to increase slightly as 
L increases. Such a feature is analogous to those of the Mott 
transitions in RHM induced by D-H binding factors.'^^'^'^''^^ 

Next, we check the order of this transition. In Fig. 3, the 
quasiparticle renormalization factor Z is shown vs U/t; Z is 
obtained using Z = ^(^p - 0) - ^(^f + 0) on the F-X seg- 
ment, where the values of n(k) at k ^ k^^ ± are estimated 
using third-order least- squares fits of the data for k < k^ and 
k > k^, respectively. There exist clear discontinuities in Z at 
U = Uc. The optimized spinon-binding parameter yu plotted 
in the inset of Fig. 3 also exhibits a large jump at ^ = Uq. 
In fact, other physical quantities show a similar discontinuous 
behavior. Thus, we may safely conclude that this transition is 
not a continuous transition but a first-order phase transition. 
The reason why the previous study could not find the cor- 
rect transition order is that the discontinuous behavior of ^(k) 



normal 




Fig. 3. (Color online) Quasiparticle renormalization factor as function of 
\U\/t for three particle densities, estimated from discontinuities of w(k) in 
Fig. 2 and other data. The inset shows the optimized spinon-binding param- 
eter as a function of \U\/t near U = Uc for the same systems as that in the 
main panel. 



manifests itself only for L > 18. The critical value \Uc\/t only 
slightly depends on n. 

Now, we consider the features of ^n in the non-Fermi- 
liquid regime |L^| > |L^cl- As shown in the inset of Fig. 3, 
the spinon-binding parameter ji approaches unity, suggesting 
that almost all up and down spinons are paired as singlets. In 
discussing gap formation, the small- |q| behavior of the charge 
(density) and spin structure factors 



A^(q) 



^(q) 



N ^ 



e-"'-"{S)S)^,) 



(21) 



(22) 



provide us with useful information. Assuming that the lowest 
excitation occurs at q = 0, the energy gap in the spin sector 
between the ground state |^o) and the first excited state |^q) 
is given by the single-mode approximation (SMA)^^ as 



<^ql^q) 

= — lim-^— ^ 



<^ql^q) 



8^-o^(q) 
where K denotes the kinetic energy, |^q) = ^ql^o), and 



(23) 



(24) 



From eq. (23), we find that A^ vanishes if 5" (q) oc ^ for ^ ^ 0, 
whereas A^ becomes finite, if 5'(q) oc q^ for ^ ^ 0. The 
charge (density) gap Aa^ can be similarly treated. 

In Fig. 4, we show 5'(q) and A/^(q) for some values of \U\lt 
near U = Uc. In the vicinity of q = 0, as \U\/t increases, 
S (q) abruptly changes its behavior from linear to quadratic 
at U = Uc, as shown in the inset of Fig. 4(b). Thus, it is 
very likely that the spin gap is generated in the non-Fermi liq- 
uid regime. In Fig. 5, we plot the spin gap estimated using 
eq. (23) for the segment of (0, 0)-(7r, 0) of S(q); the magni- 
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Fig. 4. (Color online) (a) Spin and (b) charge (density) structure factors 
along path (0, 0)-(7r, 0)-(;r, 7r)-(0, 0) for some values of \U\lt near \Uc\lt (9.0 < 
\Uc\lt < 9.05). The inset in (b) is a magnification of 5(q) near q ~ on the 
segment (0, 0)-(7r, 0). The tendency is the same for other values of n. 
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Fig. 5. (Color online) Spin gaps estimated by single-mode approximation 
for three particle densities as functions of \U\lt. The dash-dotted line is a 
visual guide for n = 0.5. The spin-gap transition occurs at \Uc\/t ~ 9. 



tude of A^ is proportional to |^| (A^ ~ 0.13|^| for n = 0.5), 
and depends on n only weakly. In contrast, the behavior of 
N(q) shown in Fig. 4(b) is almost unchanged including the 
2^F anomaly, if |6^|/^ varies, and remains linear in ^ for ^ ^ 
for |[/| > \Uc\. Thus, low-energy properties with respect to the 
charge degree of freedom are unlikely to be affected by this 
transition; regarding the charge excitation, ^n remains gap- 
less and conductivity is preserved in the spin-gapped regime. 
This may be the first realization of a conductive spin-gapped 
normal state in the variation theory. 
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Fig. 6. (Color online) (a) Doublon density as function of \U\/t for ^n and 
^sc- The dashed lines indicate the maximum D, i.e., n/2. The inset in (b) 
shows the magnification of D (^n) for n = 0.25 near the transition point, at 
which D displays a small jump, (b) Average distances from up spinon to its 
nearest down spinon for three particle densities as functions of \U\/t. 



3. 3 Picture of transition and spin-gapped state 

To deepen our understanding of the above spin-gap transi- 
tion, let us look at some other quantities. Figure 6(a) shows 
the doublon density 



D= —Yiblhi). 



(25) 



A^\U\lt increases, D increases in the Fermi-liquid state owing 
to the attractive correlation of !Pg, but it reaches almost its full 
value {nil) at U = Uq. The main panel of Fig. 6(b) shows 
the average distance from an up (down) spinon to its nearest 
down (up) spinon r^d- Here, we measure distance r by the 
stepwise (so-called Manhattan) metric. As \U\/t increases in 
a small-IL'^l/^ regime, r^d increases because the densities of up 
and down spinons decrease owing to doublon formation, and 
the binding correlation of Pq is still weak, as in the inset of 
Fig. 3. However, r^d abruptly drops when U approaches Uc, 
and converges to unity for |6^| > \Uc\, because an up spinon 
and a down spinon are tightly bound within NN sites (^ ^ 1). 
Consequently, for |6^| > \Uc\, almost all particles form onsite 
pairs, and even if a doublon resolves into spinons, they remain 
an adjacent pair and do not itinerate as isolated spinons. 

Thus, we notice that this spin-gap transition can be under- 
stood in parallel with a recently proposed picture of Mott tran- 
sitions owing to the D-H binding.^^'^^ Here, we postulate that 
antiparallel spinon pairs with a pair domain of size ^ud are cre- 
ated by the attractive correlation of Pq. We can appropriately 
define this binding length ^ud and also the minimum distance 
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Fig. 7. (Color online) Illustration of mechanism of spin-gap transition. The 
up and down arrows denote up and down spinons, respectively. Typical con- 
figurations are shown for the two phases. Although in each pair domain drawn 
with a pale circle, at least one up spinon and one down spinon must exist, ex- 
cess spinons can still move around independently, slipping out of their origi- 
nal domain, as indicated by the long dashed arrows in the Fermi-liquid phase. 
In the spin-gapped phase, a spinon cannot itinerate independently of the part- 
ner spinon out of the pair domain. 
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Fig. 8. (Color online) The binding length of up and down spinon pairs, ^ud. 
and the minimum distance between spinon pairs, ^uu, defined by eqs. (26) and 
(27), respectively, are plotted as functions of U/t for three densities. An arrow 
indicates the spin-gap transition point, the n dependence of which is small. 
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from a spinon to its nearest parallel spinon ^^u as 



^. 



ud 



^Ud + CTud, 



Suu ~ ^uu ^u 



(26) 

(27) 



where Tuu is the average distance from an up (down) spinon 
to its nearest up (down) spinon, and cr^d and o-^u are the stan- 
dard deviations of r^d and r^u, respectively. In the spin-gapped 
phase, the relation ^ud < ^uu holds, indicating that the domains 
of pairs do not usually overlap, at least, not in sequence. Con- 
sequently, almost all pairs are isolated and an up spinon and a 
down spinon are confined within ^^d, resulting in singlet pairs 
of small lengths with finite excitation gaps. In contrast, in the 
Fermi-liquid phase, ^ud becomes longer than ^uu, indicating 
that the domains of spinon pairs overlap with one another. 
Then, an up spinon in a pair can exchange a partner down 
spinon with a down spinon in an adjacent pair. As a result, 
an up spinon and a down spinon can move independently by 
exchanging their partner, as shown in long arrows in Fig. 7, 
and definite singlet pairs cannot be specified. Thus, as \U\/t 
is varied, a spin-gap transition takes place when ^ud becomes 
equivalent to ^uu, which is expected to be a monotonically in- 
creasing function of \U\/t. 

Figure 8 shows ^ud and ^^u estimated from the VMC results 
as functions of \U\/t for three particle densities. As expected 
from Fig. 6(b), ^ud abruptly drops at L^ = Uc, whereas ^ud 
monotonically increases as \U\/t increases with a small jump 
at the transition point. As a result, ^ud and ^uu intersect each 
other 3X U = Uc for any n. Thus, the scheme illustrated in 
Fig. 7 is justified to some extent. 

Finally, we discuss the itinerancy of particles. To this end, 
it is convenient to decompose the kinetic energy £'kin into two 
parts (£'kin = El +£^2), namely, the contribution of the hopping 
processes that do not (do) change the number of doublons Ei 
(£^2),^^ as shown in the lower part of Fig. 9. In the main panel 
of Fig. 9, El, E2, and £^kin are depicted as functions of \U\/t. 
For \U\ < \Ucl both Ei and E2 contribute to Eidn because iso- 
lated spinons are independently mobile, whereas in the spin- 
gapped phase. El almost vanishes (£kin ~ ^2) for any particle 
density. In this case, the independent motion of a spinon not 
accompanied by an antiparallel spinon is strongly suppressed. 
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Fig. 9. (Color online) In the main panel, the kinetic energy £'kin and its 
components E\ and E2 are depicted as functions of correlation strength for 
two particle densities. Ei (E2) is the contribution of hopping processes that do 
not (do) change D, as shown in the illustration in the lower part. The dashed 
line for E2/t is a guide proportional to -t/\U\ (\U\ > \Uc\). 



On the other hand, the contribution of the dissociation of a 
doublon into a spinon pair and their reunion, E2, remains ap- 
preciable and is proportional to -f/\U\ for large \U\/fs. This 
point is in sharp contrast to a feature of the Brinkman-Rice 
transition^"^ derived using the Gutz wilier approximation;^^ in 
this case, the motion of particles is completely prohibited for 
1^1 > I^brI = ^\E(U = 0)1,^^ so that the state becomes insu- 
lating, and a charge (Mott) gap opens. 

The above feature of £^1 =0 (\U\ > \Uc\) for any n is dis- 
tinct from that of the D-H binding state for strongly correlated 
RHM, where Ei is ^-dependent (oc 1 - ^) near half filling. ^^ 
Thus, in doped Mott insulators in RHM, doped holes or parti- 
cles play a role of carriers, and bound D-H pairs are localized 
and not involved in conduction. On the other hand, in AHM, 
the motion of a doublon, which is composed of two single 
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then decreases for |L^| > \Uco\ as ~ f/|[/|. As we will see later, 
various properties of SC actually exhibit qualitative changes 
at approximately this |6^col/^ from a BCS type to a BEC type. 
Note that normal-state properties are deeply involved in the 
crossover;^^ \Uco\lt is affected by the spin-gap transition point 
\Uc\lt in ^N, where E^ exhibits a cusp, resulting in Uco ~ Uq. 
Recall that the normal state ^n, underlying ^sc, is a Fermi 
liquid for I L'^ I < \Uc\, but ^n becomes a spin-gapped state 
in the absence of a Fermi surface, as shown in Fig. 2 for 
\U\lt > 9.05. Namely, for \U\ > \Ucl a SC transition cannot 
be interpreted by the instability of the Fermi surface against 
an attractive interaction. In this relation, IS^E means the con- 
densation energy for \U\lt ~ according to the BCS theory, 
but A£ probably deviates from the condensation energy ob- 
served experimentally for \U\ ^ |L^coL as in the case of high-Jc 
cuprates.^^ 
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Fig. 10. (Color online) In the main panel, kinetic energy and two elements 
of {"Hqe) [eqs. (28) and (29)] normalized by J = 2^1 U are plotted in the BEC 
regime as functions oi\U\lt. Illustrated in the lower part is the relation be- 
tween a single hopping process in "TY and a doublon hopping process in "Keff . 
The inset, discussed in §4.1, shows a comparison of the doublon hopping E^ 
between the normal and SC states as a function of particle density. 



particle-hopping processes as shown at the bottom of Fig. 10, 
contributes to the kinetic energy. To check this, we calculate 
the doublon hopping and diagonal terms in eq. (10), namely, 

— - — V/ 



0.15 



^j\ 



(28) 






(29) 



with / = 2r^/|[/|, in the original Hilbert space of eq. (1). In the 
main panel of Fig. 10, we compare them with the single hop- 
ping contribution (£^kin ~ ^2) for large \U\lf^. In addition to 
the large constant contribution of E perl J, the doublon hopping 
(£^d) has an appreciable magnitude, indicating the possibility 
of transport. Thus, the normal state ^n is conductive for any 
values of \U\lt and n {i^ 1).^^ 

4. Crossover of Superconducting Properties 

In §4.1, we discuss the BCS-BEC crossover in the light of 
energy gain in the SC transition and of chemical potential, 
and show that the SC transition in the BEC regime is induced 
by the kinetic-energy gain. In §4.2, we discuss quantities that 
characterize the BCS and BEC regimes. In §4.3, we roughly 
estimate the coherence length and interpair distance, thereby 
giving an intuitive picture of the crossover. 

4. 1 Energy gain and kinetic -energy -driven transition 

First, we discuss the BCS-BEC crossover from the point 
of view of the energy difference per site between the normal 
(^n) and SC (^sc) states ^E (> 0) defined in eq. (19). In 
Fig. 11(a), the \U\lt dependence of /^E/t is shown; /S.EIt in- 
creases as ~ Qxp(-t/U) corresponding to the BCS theory for 
small \U\/t's, reaches a maximum at |(7| = \Uco\ ~ 8.7r, and 
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Fig. 11. (Color online) (a) Energy differences between normal and SC 
states as functions of \U\/t for three particle densities. The maximum point is 
If^col/^ ~ 8.7 for these n's, as indicated by a dashed line. The dash-dotted lines 
for n = 0.5 are visual guides for the two limit \U\/t -^ [oc exp(-at/\U\)] 
and ^ 00 [oc /3t/\U\] with a and/3 being constants. The kinetic and interaction 
components of AE are drawn in (b) and (c), respectively. 

It is important to analyze AE into the kinetic part A^kin and 
the interaction part A£'int (AE = AE^^ -\- AEi^i). In the BCS 
theory, a SC transition is induced by lowering Ei^i at the cost 
of smaller loss in ^kin- On the other hand, it is known in a 
large- [//r regime of RHM that a SC transition occurs by re- 
ducing £^kin with a loss in £^int.^^'^^ In the latter case, the low- 
frequency sum rule of optical conductivity cri(<jj)^^ should be 
broken, namely, high-frequency excitations in cri(oj) arise, 
because the sum of cri(oj) is proportional to -£^kin^^ on a 
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square lattice. In Figs. 11(b) and 11(c), we show A£'kin and 
A£'int, respectively, for AHM. For \U\ < |6^col, ^int (^kin) has a 
gain (loss) by the SC transition in accordance with the BCS 
theory. Meanwhile, for \U\ > \Ucol the situation is reverse; 
the SC transition is driven by a gain in kinetic energy. Corre- 
spondingly, the hopping of doublons (carriers) \Ej)\ becomes 
more enhanced in ^sc than in ^n in the BEC regime, as 
shown in the inset of Fig. 10. Kinetic-energy-driven (SC or 
magnetic) transitions may be rather general in strongly cor- 
related systems.^^'^^'^^ As mentioned previously, this reversal 
of driven force will be experimentally found if cri (oj) is accu- 
rately measured in cold-atom systems. In fact, similar results 
were reached for AHM in infinite dimensions by DMFT.^^'^"^ 
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Fig. 12. (Color online) Optimized variational parameters, (a) Gutz wilier 
(onsite) factor g, and (b) NN antiparallel-spinon factor ju for ^n and ^sc are 
shown as functions of \U\/t. The symbols are common in (a) and (b). The 
arrows indicate the spin-gap transition points in ^n- 

Let us consider this SC transition in the light of the vari- 
ational parameters. At the level of one-body wave function, 
Obcs improves the energy over Op by creating onsite Cooper 
pairs through the pair potential Ap (Fig. 21), in accordance 
with the BCS theory. Therefore, in the BCS regime, the num- 
ber of doublons is expected to be more in the SC state than 
in the normal state. This is actually shown in Fig. 6(a). Since 
the increase in D hinders the motion of particles, the kinetic 
energy is suppressed in the SC state, as in Fig. 11(b). As 
|L^| increases, however, D of ^n increases more rapidly and 
surpasses that of ^sc at l^col- This reversal is brought about 
mainly by the correlation factor P. In Fig. 12, we compare the 
optimized parameters in P between the normal and SC states. 
The onsite attractive factor g is certainly larger in ^n, espe- 
cially near U = Uco- The antiparallel-spinon binding factor// 
mainly works for the suppression of overgrown Dby g in the 
BEC regime in order to gain £^kin. 

In a one-body framework, when the chemical potential (f ) 
is situated in an energy band, low-energy excitation in the SC 
phase is described by Bogoliubov's quasiparticles, whereas 
when f becomes lower than the band bottom si^, the statis- 
tics of the system becomes bosonic. Thus, the BCS-BEC 
crossover point is roughly estimated using f = ^p. We esti- 
mate f for ^sc from f = dE^ddn (strictly finite diff'erences). 
Within statistical error, £^sc is almost a linear function of n for 
0.5 < ^ < 1.0, so that the f obtained in this range becomes 
independent of ^. In Fig. 13, we plot the thus-estimated f as a 
function of \U\lt with the value at half filling i.e., f = Ujlt. f 
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Fig. 13. (Color online) The chemical potential estimated from Ejt at 
higher particle densities (0.5 < n < 1.0) is shown by dots as a function of 
\U\lt. The dash-dotted line denotes the value at half filling: ( = U/2. The 
dashed line indicates the band bottom sl (-40- 



reaches the band bottom si^ = -4t at \U\/t ~ 7.9. The behav- 
ior of f here is consistent with those obtained by DMFT.^^'^^ 
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Fig. 14. (Color online) Evolution of momentum distribution function for 
^sc as \U\/t varies. The dashed line indicates \U\/t = oo. 

Finally, let us look at the momentum distribution function 
for ^sc (Fig- 14). As \U\/t increases, a step-function form of 
\U\/t = changes to a BCS type (v^) for small \U\/t and is 
then smoothly modified through \Uco\/t to a constant in the 
BEC limit (\U\/t = oo). Such evolution of n(k) has already 
been observed in an experiment on an ultracold gas in a trap;^^ 
we hope for similar experiments on optical lattices. 

4.2 Pair correlation function and helicity modulus 

As mentioned in §1, a previous study of Toschi et al}^ us- 
ing DMFT argued that appropriate quantities that trace the 
strength of SC {Tc) in the BCS and BEC regimes are the gap 
Asc ~ ^AA) ^^^ th^ superfluid stiff'ness /)§, respectively. Asc 
indicates the cost of creating a Cooper pair, while D^, char- 
acterizes the cost of realizeing phase coherence. In this sub- 
section, we start with the quantities corresponding to Asc and 

A. 

As an appropriate index of off'-diagonal-long-range order 
(ODLRO) in the present scheme, we consider a SC correlation 
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Fig. 15. (Color online) The onsite superconducting correlation function is 
drawn along a path of r^, (0,0)-(L/2, 0)-(L/2, L/2)-(0,0), for various values 
of \U\lt. For other values of n and L, the behavior is basically the same. 
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Fig. 17. (Color online) Onsite superconducting correlation function of ^sc 
plotted as a function of\U\lt for three particle densities and L = 20. The dash- 
dotted line is a visual guide of oc exp(-t/\U\) for n = 0.75. 
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Fig. 16. (Color online) System size dependence of Poo for three values of 
\U\/t atn = 0.5. In the cases of L not satisfying the closed-shell condition at 
n = 0.5, we adopt the average of n = 0.5 ± S with the smallest S satisfying 
the condition. As visual guides, we show the first-order least-squares fits for 
L ^ oo with dashed lines. Solid (open) symbols represent the data of ^sc 
(^n)- ^n is bound to vanish for L ^ oo. 



function of onsite pairing, ^^ defined as 

P(r,) = l^<M/,^,,). 



(30) 



The magnitude of ODLRO is given by the long-distance value 
of P(r^), i.e.. Poo = limir^i^co Pir^) ~ A^^^. In the present VMC 
calculations with finite systems, we must check the r^ depen- 
dence of P(r^). In Fig. 15, we plot P(r^) for various \U\/fs 
along a typical path on the lattice. Since P(r) is substantially 
constant for |r| > 2 for any value of \U\/t, consistently with a 
QMC study,^^ it is appropriate to put P(r) with the most dis- 
tant r = (L/2, L/2) at Poo, and check its system-size depen- 
dence. In Fig. 16, we plot the thus-estimated Poo as a function 
of 1/L for three values of \U\/t. In the SC state, the system- 
size dependence of Poo is very weak for any \U\/t, and fitted 
well by a first-order least-squares method. Thus, we may dis- 
cuss Poo with a finite but large L. Figure 17 shows the \U\/t 
dependence of Poo for L = 20.IntheBCS regime (|^|<| ^co I), 
Poo increases as ~ exp(-t/U) as \U\/t increases, as predicted 
by the BCS theory. In this regime, the optimized Ap is consid- 
ered to approximate the gap, as will be argued; we confirmed 



using the data shown in Fig. 21 that the relation Poo ~ Ap 
holds. On the other hand, in the BEC regime. Poo tends to con- 
verge at a finite value as \U\/t increases, in accordance with 
the result of Asc obtained using DMFT.^^'^"^ Thus, the behav- 
ior of Poo or Asc in this regime does not coincide with the 
behavior of Tc, ~ t^/\Ul which was naturally expected^ and 
actually obtained by DMFT.^^'^^'^^ 

Incidentally, a similar SC correlation function with the NN 
dx2_y2-w2iYe pairing PJ has been calculated for RHM in 2D by 
VMC with the same class of trial functions. ^^ In the strongly 
correlated regime (typically U/t = 12), where the cuprates 
are considered to be properly described, PJ behaves as the 
so-called dome shape as a function of doping rate S (= l-n). 
This dome shape closely agrees with the S dependence of T^ 
experimentally observed for the cuprates. If the framework of 
BSC-BEC crossover as a function of S is applicable to this 
case, PJ decreases and scales with Tc as the parameter ap- 
proaches the BEC limit (^ ^ 0). Consequently, the behavior 
of P J in RHM does not fully correspond to that of Poo in this 
study. 




Fig. 18. (Color online) Helicity modulus as a function of \U\/t for three 
particle densities. The dash-dotted line indicates a visual guide of oc -t/\U\ 
for ^ = 0.75. 

Now, we turn to the helicity modulus ps, which is related to 
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superfluid stiffness D^ with ps = D^IAne^. In the SC state, D^ 
is equivalent to the strength of the delta- function component 
at 6l) = in the optical conductivity cri {cSf^ and represents 
the superfluid weight. We calculate ps as the coefficient of the 
quadratic term in the increment in energy when the phase of 
order parameter A^ is twisted by q as Ay = |A|^'*i'^^:^^ 

E(q)-£(0) = 2psq'+O(q^ (31) 

following the prescription of Denteneer et alP^ In Fig. 18, the 
\U\lt dependence of the ps thus obtained is plotted for three 
particle densities. The resultant ps here is almost independent 
of the q used (|q| ~ 0.1), and the system size dependence is 
negligible within the symbols between L = 12 and 20. Re- 
gardless of n, Ps is a monotonically decreasing function of 
\U\lt. The behavior in the BCS regime is distinct from that of 
Tc ~ exp(-^/|^|), butps scales with T^ ~ f- l\U\ in the BEC 
regime, indicating the strength of SC in the BEC regime is 
determined not by the cost of creating a pair but by the cost 
of realizing phase coherence. The present result of ps is con- 
sistent with the previous results obtained by a VMC method 
with ^G<^BCs7^ QMC,20 and DMFT.23'24,26 




Fig. 19. (Color online) Spin gap of SC state estimated similarly to that in 
Fig. 5 by single-mode approximation. The dash-dotted straight line indicates 
an extrapolation from large- |f/|/f values for n = 0.5. 



In the remainder of this subsection, we discuss some quan- 
tities related to P(r) and ps. First, we take up the small- |q| 
behavior of spin and charge structure factors, eqs. (21) and 
(22). Like for the normal state, A^(q) oc |q| for |q| -^ for 
any \U\/t, showing that ^sc is conductive in particle density. 
On the other hand, S(q) oc |qp for any \U\/t (> 0) in contrast 
to the case of ^n, indicating that a spin gap opens owing to 
pair formation. We estimate the spin gap A^ for ^sc using 
SMA [eq. (23)], and show the \U\/t dependence in Fig. 19. 
Although we do not display the data for small \U\/fs owing 
to the relatively large errors due to the use of a finite system 
(L = 20), A^ seems to be proportional to Qxp(-at/\U\) for 
small \U\/fs. On the other hand, A^ is proportional to \U\/t, 
and has a magnitude similar to that of ^n- ^s is almost inde- 
pendent of ^. Thus, A^ is a quantity that scales with Tc in the 
BCS regime. 

Next, let us consider the condensate fraction po. We may 
regard b^. as a creation operator of a hard-core spinless boson 




Fig. 20. (Color online) Condensate fractions of hard core bosons (or dou- 
blons) for three particle densities as functions of \U\/t. The dash-dotted line 
is a visual guide of oc exp(-t/U) for n=0.26. 



at the site /'; b^. satisfies the Bose commutation relation ex- 
cept for the same site. Following Bose systems,"^"^ we define a 
quantity corresponding to the condensate fraction or the k = 
element of momentum distribution function hdO^) for b : 

p, = ^nom=^Yj<^]bj,,). (32) 

In Fig. 20, the \U\/t dependence of po is depicted for three 
particle densities. The behavior of po for small \U\/fs is 
Po ~ exp(-t/\U\) and has a meaning similar to Poo. In the BEC 
regime (\U\ > \Uco\), Po is almost constant, indicating that a 
picture of hard-core bosons is justified in the entire regime of 
BEC. The suppression of po with increasing n is primarily be- 
cause a high density enhances the effect of onsite interaction. 




n, L 

-*— 0.75,20 

0.63, 20 

-•—0.5, 20 

o 0.37,20 



—0.26, 20 
-^—0.13,20 

^G^BCS 

0.5, 12 



'" mf, "• 



Fig. 21. (Color online) Optimized gap parameter Ap for several «'s as func- 
tion of \U\/t. The dash-dotted line is a visual guide of oc Qxp(-t/\U\). The inset 
shows a comparison of the optimized values of Ap between IPqPg^bcs and 
'^G<I>BCsfor^ = 0.5. 

Finally, we consider the pairing gap parameter Ap given in 
eq. (16). In Fig. 21, we show the \U\/t dependence of the op- 
timized Ap. For \U\ < \U\co, it is natural that Ap represents the 
SC gap of oc exp(-t/\U\), as expected from the BCS theory. 
In the BCS theory, Asc should continue to linearly increase 
like A^ in Fig. 19. However, the Ap of ^sc exhibits a peak 
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at L^ ~ Uco, and then decreases as \U\lt increases, similarly 
to A^" [Fig. 11(a)]. In the BEC regime, it is probable that Ap 
obeys oc 2fi/\U\ = /, because / is the sole energy scale for 
\U\/t -^ oo, according to eq. (10). It seems that AE and Ap 
scale Tc in AHM. 

In the inset of Fig. 21, we compare the optimized Ap's be- 
tween ^sc and !Pg^bcs- Although the two Ap's behave sim- 
ilarly in the BCS regime, the Ap of !Pg^bcs monotonically 
increases unlike the Ap of ^sc in the BEC regime. It fol- 
lows that the binding correlation between antiparallel spinons 
is also significant for ^sc, especially in the BEC regime. 

In variational theories with J;^2_3;2-wave SC states for 
cuprates, the J;c2_^2-wave gap parameter A^, corresponding 
to Ap here, is considered to represent a singlet-pairing gap 
(not necessarily SC gap).^"^'^^ The optimized A^ monotoni- 
cally increases as the doping rate, the relevant parameter of 
the crossover, approaches the BEC limit (^ ^ 0), in contrast 
to TcP~^^ The behavior of Ap here is distinct from that of 
A^. Again, we should be careful to consider the cuprate in the 
point of view of the BCS-BEC crossover. 



iMj^ 




\U\/t 



Fig. 22. (Color online) Two estimates (a/Ap and |ud) of the pair size ^pair, 
corresponding to the coherence length, and the average minimum distance 
between Cooper pairs |uu are compared as functions of \U\/t. The dashed 
line indicates \U\/t where ^pair intersects ^uu- 



4.3 Coherence length and intuitive picture 

The BCS-BEC crossover has been typically explained by 
whether or not a domain of a Cooper pair overlaps with a do- 
main of another pair, as in Fig. 7. To discuss this more quan- 
titatively, we need to estimate a pair size ^pak corresponding 
to the coherence length and a distance between Cooper pairs 
luu. As for ^pair, it is reasonable to refer to the BCS expression 
of Pippard's coherence length. 



^0 = 



(33) 



^lAscI 

in the BCS regime. In the present study, vp is a constant for 
any \U\lt, because the renormalization of ^p by \U\lt is not 
introduced. Thus, we assume ^pak = Q^/Ap, where a is a con- 
stant determined so that ^pair can be smoothly connected to the 
form on the BEC side. In the BEC regime, eq. (33) does not 
work, because vp cannot be defined. Thus, following eq. (26), 
we naively assume ^pair = ^ud = ^ud + <^ud. where r^d and a^^ 
denote the average distance between a spin (not necessarily 
of a spinon) and its nearest antiparallel spin and the standard 
deviation of rud, respectively. Note that Tud (o^ud) is diff'erent 
from Tud (cTud) in eq. (26) in that spins that constitute doublons 
are taken into account. Similarly, following eq. (27), we esti- 
mate an interpair distance as the average minimum distance 
between a spin and its nearest parallel spin, |uu = f^^ - d-^^. 

In Fig. 22, the \U\lt dependences of ^pak and |uu thus es- 
timated are compared for n - 0.5. The pair size ^pak is a 
monotonically decreasing function of |6^|/^, whereas the inter- 
pair distance ^uu is almost independent of 1 1/|/^. Consequently, 
^pak crosses luu at 1^1 = \l]^\ ~ 6.2^ at this particle concentra- 
tion. Thus, for \U\ < \U^\, Cooper pairs penetrate each other 
(^pak > luu) as in the Fermi liquid phase in Fig. 7. On the other 
hand, for IL'^I > \U^\, a pair becomes almost a point hard-core 
boson (^pak ~ 0), and is isolated from other pairs (^pak < luu)- 
Similarly, we estimated \lJ^\lt for other values of ^, and found 
\U^\ ~ 6^ irrespective of particle densities. ^^ Since the above 
estimation of U^ is rather broad, we consider that U^ should 
be identical to L^co- 

Finally, we point out the diff'erence in pairing manner be- 
tween the spin-gap transition at U^ in the normal state (§3.3) 



and the crossover at ~ U^q in the SC state. Bound spinon 
pairs in the normal state {\l]\ > \Uc\) dissociate into indepen- 
dent spinons immediately when the pair domains overlap with 
each other dXU = 11^. On the other hand. Cooper pairs in the 
SC state remain paired, even if pair domains come to consid- 
erably overlap for \U\ ^ \Uco\ in the BCS regime; there is 
no critical change at t/ = t/co- Consequently, a phase transi- 
tion (crossover) arises and a spin gap closes (survives) on the 
weakly correlated side in the normal (SC) state. The stability 
of Cooper pairs against mutual overlap is a current topic. ^^ 

5. Conclusions 

Using a variational Monte Carlo (VMC) method, we stud- 
ied the features of a spin-gap transition in a normal state and 
of the BCS-BEC crossover in a superconducting (SC) state 
in the attractive Hubbard model (AHM) on the square lattice. 
We summarize the main results below. 

(1) In the normal state, we revealed that, unlike the sim- 
ple Gutz wilier wave function (GWF), a wave function with 
an antiparallel- spinon binding correlation Pq [eqs. (12) and 
(13)] undergoes a first-order transition from a Fermi liquid to 
a spin-gapped phase at \Uc\lt ~ 9. In the spin-gapped phase, 
particle density current can flow through the hopping of dou- 
blons. The pseudogap phase above T^ for |t/| > |6^cl may be 
deduced from the properties of this wave function. The mech- 
anism of this spin-gap transition is understood to be simi- 
lar to that of a Mott transition in a repulsive Hubbard model 
(RHM) induced by a doublon-holon binding correlation."^^' "^"^ 
We would also like to realize a variational normal state that is 
spin-gapped and conductive in RHM. 

(2) We first appHed VMC to the SC state of AHM, and con- 
firmed that, as \U\lt increases, the mechanism of supercon- 
ductivity undergoes a crossover at approximately |^col ~ l^d 
from an BCS type to a Bose-Einstein condensation (BEC) 
type. Pq is again needed to suppress the gap, which is greatly 
overestimated in GWF for \U\ ^ |L^col- In the weak-correlation 
regime (\U\ < \Uco\), the strength of SC (T^c) is scaled with 
quantities related to the SC gap as ~ exp(-r/|6^|), as expected 
from the BCS theory. For \U\ > \UcoV the superfluid stifl'- 
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ness, which is related to the cost of phase coherence, scales 
with Tc as fi/\U\. Such typical features of this crossover are 
captured by the energy gain in the SC transition AE in the 
whole range of \U\/t. In the BEC regime, the SC transition 
is induced by a gain in kinetic energy; this aspect is in con- 
trast to the BCS theory, but is in accord with the magnetic and 
SC transitions in strongly correlated RHM."^^'^^ Most features 
are consistent with the framework of BCS-BEC crossover that 
previous studies provided. 

(3) The physics of a spin-gap transition in the normal state 
and the BCS-BEC crossover in the SC state are explained in 
a semiquantitative manner by a simple idea of the competi- 
tion between the pair size ^ud and the interpair distance ^uu, as 
shown in Fig. 7. This idea is equivalent to that of Mott transi- 
tions in RHM,'^^'^'^ in which a doublon-holonpair corresponds 
to the singlet pair here. 

(4) In connection with that observed high-Tc cuprates, the 
\U\/t dependence of the pair correlation function Poo and 
the gap parameter Ap studied here qualitatively differ from 
the doping rate (S) dependence of the corresponding quan- 
tities (PJ and Ad) in the strongly correlated RHM, when 
the relevant parameters (\U\/t and S) are in the respective 
BEC regimes. Furthermore, in a strongly correlated RHM, the 
dx2_y2-w8iYe SC transition is always kinetic-energy-driven, re- 
gardless of nP We will address this subject more carefully in 
upcoming publications. 
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